Kinetic theory of cavity cooling and self-organisation of a cold gas 
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We study spatial self-organisation and dynamical phase-space compression of a dilute cold gas of 
laser-illuminated polarisable particles in an optical resonator. Deriving a non-linear Fokker-Planck 
equation for the particles' phase-space density allows us to treat arbitrarily large ensembles in 
the far-detuning limit and explicitly calculate friction forces, momentum diffusion and steady-state 
temperatures. In addition, we calculate the self-organisation threshold in a self-consistent analytic 
form. For a homogeneous ensemble below threshold the cooling rate for fixed laser power is largely 
independent of the particle number. Cooling leads to a g-Gaussian velocity distribution with a 
steady-state temperature determined by the cavity linewidth. Numerical simulations using large 
ensembles of particles confirm the analytical threshold condition for the appearance of an ordered 
state, where the particles are trapped in a periodic pattern and can be cooled to temperatures close 
to a single vibrational excitation. 

PACS numbers: 37.30.-|-i, 37.10.-x, 51.10.+y 



I. INTRODUCTION 

A dilute cold gas of polarisable particles can be manip- 
ulated in a controlled way using the light forces induced 
by a sufficiently strong laser far off any internal optical 
resonance [1]. For a free-space laser setup this force gen- 
erates a conservative optical potential for the particles 
with a depth proportional to the local light intensity. As 
the forces are generated via photon redistribution among 
different spatial directions, the particles in turn alter the 
field distribution and act essentially as a spatially varying 
refractive index. While this backaction can safely be ig- 
nored in standard optical traps [2] , it was shown to have 
a significant effect if the light fields are confined within 
an optical resonator enhancing the effective particle- light 
interaction [3]. 

For transverse illumination a threshold pump inten- 
sity where this coupled particle-field dynamics can lead 
to spatial selfordering of the particles into a regular 
pattern — resembling very closely a phase transition — was 
theoretically predicted and experimentally confirmed [4- 
6]. Due to cavity losses this dynamics is dissipative and 
thus can constitute a new cooling mechanism for a very 
general class of polarisable objects [7, 8] . Extensive sim- 
ulations using fairly large numbers of particles with large 
detunings predict that already with current molecular 
sources and cavity technology a useful phase-space com- 
pression could be achieved [9]. However, particle-based 
simulations cannot be applied to sufficiently large particle 
numbers and laser powers for the whole parameter range 
of interest. As an alternative, a Vlasov-equation-based 
approach [10] for the particles' phase-space distribution 
provides a description for arbitrarily large ensembles, but 
excludes correlations important on longer time scales. In 
this work we adopt methods from plasma kinetic theory 



to generalise the Vlasov model to include correlations, 
leading to a non-linear Fokker-Planck equation for the 
statistically averaged phase-space distribution, which in- 
cludes friction and diffusion and allows to predict cooling 
time scales and the unique steady-state distribution. 



II. SEMICLASSICAL EQUATIONS OF MOTION 

We consider N polarisable particles moving along the 
axis of a lossy standing- wave resonator assuming strong 
transversal confinement. The particles are off- resonantly 
illuminated by a transverse standing-wave pump laser 
and scatter light, whose phase is determined by the 
particle positions, into the cavity. The quantum mas- 
ter equation describing this system can be transformed 
into a partial differential equation for the Wigner func- 
tion [11]. Truncating this equation at second-order 
derivatives (semiclassical limit) yields a Fokker-Planck 
equation, which for positive Wigner functions — excluding 
all non-classical states [12] — is equivalent to the Ito 
stochastic differential equation (SDE) system 



dxj = Vj dt 

1 dU(xi,a) , 
dvj = 1 ^' ' dt 



dxj 



(la) 
(lb) 



da = 



N 



i Ac -C/o^sin2(fca;j) 



idt~ 



N J— 

- ir^^sinikxj) dt+J^ {dWi +idW2), (Ic) 

with the single-particle potential 

U{x, a) = hUo\a\'^ sin^(fcx) -|- hri{a + a*)sm{kx). (2) 
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Each particle has the mass m; Xj and Vj denote the 
centre-of-mass position and velocity of the jth particle. 
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FIG. 1: (Colour on-line) An ensemble of particles is illumi- 
nated by a transverse standing-wave laser and scatters light 
into the resonator (effective pump strength rj). Above a 
threshold pump intensity the particles selforganise in a pe- 
riodic pattern. Cavity losses are characterised by the decay 
rate k. 



respectively. Uq is the light shift per photon and k the 
cavity decay rate; the associated input noise is taken into 
account by the two Wiener processes dM^i and d W2 . Here 
we have neglected momentum diffusion caused by spon- 
taneous emission, which is valid for large ensembles and 
large detunings [8]. The transverse laser standing wave 
gives an effective position-dependent pump of magnitude 
77 of the field mode a with wave number k. This laser 
is detuned by Ac = Wp — Wc from the bare cavity res- 
onance frequency. We will focus on the weak-coupling 
limit iV|t/o| <C |Ac| ~ k throughout this paper. Refer to 
fig. 1 for a schematic view of the system. 

We used the scheme proposed in [13, 14] for the di- 
rect numerical integration of the SDE system (1). How- 
ever, for analytical predictions and the description of very 
large ensembles, a continuous phase-space description as 
described below proves far more suitable. 



III. CONTINUOUS DESCRIPTION, 
INSTABILITY THRESHOLD 



The semiclassical SDEs 
Klimontovich equation [15] 



(1) are equivalent to the 



dfK ^ ^ dfK 1 dUdfK 
dt dx m dx dv 







(3) 



together with an evolution equation for the field mode 
a obtained replacing the sums in eq. (Ic) by the inte- 
grals A'' JJ •fK{x,v,t)dxdv. /k(x, is the so-called 
Klimontovich or "exact" distribution function with ini- 
tial condition 

1 ^ 

/k(x, t;, 0) - - ^ <S(a; - x,{0))S{v - v,iO)) . (4) 

As this function is highly irregular, the above reformula- 
tion has no computational merit by itself, but provides 



an ideal starting point for a statistical treatment. To this 
end we decompose every quantity (/k, a and U) into its 
smooth mean value and fluctuations. 



/k(x, V, t) =: (/(x, w, t)) + Sf{x, V, t), 



(5) 



with {6f{x,v,t)) = 0. The statistical average (•) is over 
an ensemble of similar initial conditions {{xj(0),Vj{0))} 
and a{0), as well as over the realisations of the white 
noise process mimicking the input noise for the cavity 
field. 

For the smooth ensemble-averaged Klimontovich dis- 
tribution, called one-particle distribution function, we 
then find 



d{f) , d{f) 1 d{u)d{f) 



dt 



dx 



m 



dx dv 



/d5UdSf\ 



Neglecting all correlations in eq. (6) leads to the Vlasov 
equation, which becomes exact in the limit N ^ 00 [10]. 

The Vlasov equation together with the equation for 
(a) possesses an infinite number of possible spatially ho- 
mogeneous stationary solutions with zero cavity field, of 
which, however, not all are necessarily stable against 
perturbations. Indeed, for any (dimensionless) sym- 
metric velocity distribution g{v/vT) := Lwt (/(w)) and 
(5 Ac - NUo/2 < we find, that if 



vp 



-2C 



(7) 



where vp denotes the Cauchy principal value, perturba- 
tions trigger a self-organisation process. Here we have 
defined the thermal velocity 2kBT/m; L is the 

cavity length. For a Gaussian distribution the integral 
evaluates to one. The relation (7) has been derived by 
methods presented in [10]. There it was shown that the 
threshold can be computed from the zeros of the disper- 
sion relation (Re(s) > 0) 



D{.s) = {s + Kf + 6^ 
.NLrf 



ihkS- 



2m 



d. if) d, if) 



ikv 



ikv 



dv. (8) 



IV. KINETIC EQUATION FOR THE VELOCITY 
DISTRIBUTION IN THE NON-ORGANISED 
PHASE 

Let us now return to eq. (6), which for weak spatial 
inhomogeneity — as expected below the instability thresh- 
old (7) — approximately reads 



d{f) /d6U d6f\ 
\ dx dv I 



dt 



(9) 



where the overbar denotes the spatial average. The right- 
hand-side correlation function can be computed using es- 
tablished methods from plasma physics as in [15] when 
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the fluctuations evolve on a much faster time scale than 
the average values, which are considered "frozen". This is 
justified as long as the system remains far from instabil- 
ity. After some lengthy calculations, which we omit here 
for the sake of compactness, we obtain (for symmetric 
distributions) the non-linear Fokker-Planck equation 



A[{F)] (F) 



|;U(^)]|;(^>) (10) 



for the velocity distribution {F{v,t)) 
the coefficients 



L{f{v,t)), with 



A[{F)] := 
B[{F)] := 



2hk5mf' kv 

m |D(ifcw)|2 
h^Prf'K + (52 + k'^v'^ 



2inn? 



\D{ikv)\'^ 



(11a) 
(lib) 



This equation, which describes the sub-threshold gas 
dynamics, is one of the central analytical results of 
this work. Note that the coefficients functionally de- 
pend on {F) through the dispersion relation D{\ijj) :— 
lim£j_o-D(e + ii^). A and B represent the deterministic 
part of the equations and the field noise, respectively. 
All cavity-mediated long-range particle interactions are 
encoded in the dispersion relation. Note that very far 
below threshold the full dispersion relation (8) reduces 
to D{ikv) ~ {ikv + k)^ + 5"^, which corresponds to the N 
independent-particles case. 



EQUILIBRIUM DISTRIBUTION AND 
TEMPERATURE 



Normalisable steady-state solutions of eq. (10) exist 
only for negative detuning S < 0, where light scatter- 
ing is accompanied by kinetic-energy extraction from 
the motion. For positive detuning the particles are 
heated. Interestingly, we obtain the non-thermal q- 
Gaussian velocity distribution function (also known as 
Tsallis distribution and closely related to Student's t- 
distribution) [16, 17] 



{F{v))^ i-(i-g) 



\ 1-? 



2kBT, 



(12) 



with q := 1 + and the recoil frequency wr :— 

hk^ /2m. We have defined an effective "temperature" pa- 
rameter 



~2' 



(13) 



The latter minimum value of T appears for 5 = —n. The 
parameter |(5|/wr determines the shape of the distribu- 
tion and gives rise to further restrictions. Normalisable 
solutions exist for \5\ > ci;p;/2, the second moment (ki- 
netic energy) only for \5\ > 3ajR/2. The case \S\ = ojr cor- 
responds to a Lorentzian distribution. For I^I/cjr oo. 
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FIG. 2: (Colour on-line) Normalised g-Gaussian velocity dis- 
tribution for S = ~2.5ijjR obtained by integrating the semi- 
classical equations (1) up to ujjit = 5N (solid line) compared 
to the theoretically predicted g-Gaussian (12) (circles). The 
Gaussian (dashed line) is plotted for the temperature (13). 
The kinetic temperature differs considerably from eq. (13), 
ftBTkin = 2.5fcBr. Parameters: iV = 5000, NUo = -O.Iljr, 
\/iV»7 = ISOOtJR, K = 100a;R and Ac = -2.55a;R. Ensemble 
average over 25 initial conditions and 10 realisations of the 
white noise process. 



i. e. g — > 1, the distribution (12) converges to a Gaussian 
distribution with kinetic temperature k^Tj^in '■— ralv'^^ 
given by the parameter k^T defined in eq. (13), which 
justifies the choice of "temperature" above. Refer to 
fig. 2 for an example of prominent g-Gaussian behaviour 
{q = 1.4). Tsallis distributions have already been ob- 
served experimentally in dissipative optical lattices [18]. 

Of course the distribution (12) can be unstable and in 
fact is, if 



2\S\ 



s=- 



■3-q 



3-g 



(14) 



which can be seen by inserting the g-Gaussian into (7). 
For a Gaussian and S = —k this criterion may be refor- 
mulated as 



N\Uo\Vopt > hK^ 



(15) 



where Vo-pt is the optical potential depth created by the 
pump laser. If the condition (14) is satisfied, we conclude 
that there exists no spatially homogeneous steady state 
at all and consequently we expect the system to organ- 
ise. The possibility of attaining such an inhomogeneous 
equilibrium even though the uniform distribution (12) is 
stable cannot be ruled out mathematically but seems un- 
likely on physical grounds. These expectations are con- 
firmed by all our numerical simulations of eqs. (1). 

Based on these considerations we predict the occur- 
rence of dissipation-induced self-organisation if the self- 
consistent relation (14) is fulfilled. That is, any initially 
stable, unorganised distribution will loose kinetic energy 
(cavity cooling) and eventually selforganise. Contrary 
to the self-organisation process of an initially unstable 
state, which is abrupt and accompanied by strong heat- 
ing, dissipation-induced self-organisation is characterised 
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by a much slower buildup of the photon number and 
monotonous cooling. 



VI. SELFORGANISED PHASE 

Above threshold the mathematics becomes more com- 
plex due to the inhomogeneous spatial distribution. How- 
ever, using action-angle variables [19, 20], we can still 
derive a Fokker-Planck equation similar to eq. (10) [21]. 
In the limit of deep trapping, where a harmonic approx- 
imation for the potential becomes valid, the steady-state 
solution for the strongly organised phase is a thermal dis- 
tribution with a kinetic temperature depending on the 
linewidth k and the trap frequency wq. 



The trap frequency is ujq 
approximated by 



+ 52 + 4a;2 



(16) 



4:r]ujY{ (I Reofool) and can be 



N 



(17) 



in the far-detuned regime where \S\ » wr; iJc is the self- 
consistent critical value defined in eq. (14). As the tem- 
perature depends explicitly on the laser power, higher 
pump strengths result in deeper trapped ensembles with 
increased kinetic energy. Note that this system has the 
interesting property that the more particles we add, the 
deeper the optical potential gets as is the case for self- 
gravitating systems [22]. 
The order parameter 



e := lim 



sin{kx) (/(x, V, t)) dxdv (18) 



is an adequate measure of particle localisation in equilib- 
rium as it is zero for a completely homogeneous distribu- 
tion and plus/minus one for the perfectly selforganised 
phase (i.e. J-peaks), the sign depends on whether the 
odd- or even wells are populated [6]. In fig 3 we de- 
pict its behaviour as a function of the pump strength. 
Initially, the particles were spatially homogeneously dis- 
tributed with a (stable) Gaussian velocity distribution; 
self-organisation sets in because of the cavity cooling 
effect. Hence the branching point is given by the self- 
consistent threshold (14). 

For the strongly organised phase, where the distribu- 
tion function is a thermal state with temperature (16), 
the order parameter reads (in harmonic oscillator approx- 
imation) 



e = ± 1 



(19) 



Its maximum value is limited by the detuning and the 
recoil frequency, & — >■ ± (1 — wr/Ii?!) for ry — 00. This 
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FIG. 3: (Colour on-line) Order parameter (18) as function of 
the pump strength obtained from eqs. (1) in the long-time 
limit {idnt = 50A*'). The red dashed line is given by eq. (19), 
the black dashed line around the critical point corresponds to 
eq. (21). Parameters: A*" = 1000, NUo = -cjr, k = IOOojr, 
Ac = NU0/2-K (i.e. 5 = -k) and k^To = llQEn- Ensemble 
average of 20 (away of the critical point) and 60 (around the 
critical point) noise trajectories, respectively, for one initial 
condition. 



corresponds to a Gaussian spatial distribution with width 
fc2(Ax)2 = 2wr/|(5|. Let us also briefly investigate the 
opposite limit of pump strengths slightly above threshold. 
Solving the self-consistency equation 



-Nt^\5\ 
;2 -f ,52 



e 



(20) 



perturbatively around the critical point for a thermal 
state yields 



6: 



±2W^ 



(21) 



and thus a critical exponent of 1/2, as already predicted 
in [6]. 

The self-consistent phase diagram including cavity 
cooling is sketched in fig. 4. For small [^I/wr the 
ensemble will always remain weakly organised (i. e. a 
small spatial modulation on top of a homogeneous back- 
ground) above threshold as the necessary prerequisites 
for a strongly organised equilibrium — and hence the va- 
lidity of eq. (16) — cannot be fulfilled. 



VII. COOLING TIME 

Let us take a closer look at the cooling time r, i. e. 
the time characteristic for the kinetic-energy equilibri- 
sation. First we treat the case of a fixed ensemble size 
and variable pump strength. The drift term A + dyB in 
the non-linear Fokker-Planck equation (10) — scaling as 
~ if/uJK — might suggest the conclusion that the larger 
the pump strength the shorter the cooling time. However, 
numerical simulations (cf. fig. 5) prove this expectation to 
be somewhat misleading. The reason therefor is the onset 
of self-organisation which occurs as soon as the momen- 
tary distribution becomes unstable according to eq. (7). 
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FIG. 4: (Colour on-line) Schematic view of the phase diagram 
in the weak-coupling limit {N\Uo\ <C ac) for k = IOOojr. Equi- 
librium solutions exist only for S < —iJB./2, the Lorentzian 
corresponds to the case \5\ — lj-r,. For large values of the 
detuning \5\, strongly organised equilibria exist already for 
pump strengths slightly above the critical value, cf. also 
fig. (3). 



Particle trapping is a hindrance to optimal cooling in a 
twofold way. Firstly, with its appearance kinetic energy 
is dissipated at a significantly lower rate because a part 
of the laser power is utilised for the buildup of the po- 
tential and consequently is not available for cooling. Sec- 
ondly, the lowest achievable temperature increases with 
the pump strength, cf. eq. (16). The situation is worst for 
initially unstable ensembles since the particles are heated 
before being cooled. Hence it is a plausible requirement 
that self-organisation has to be avoided to realise the 
optimal cooling time. Accordingly, as a rule of thumb 
we may state that the latter is achieved for a laser power 
which renders the desired gas temperature critical. Hence 
the optimal cooling time is estimated from eq. (10) to be 



'opt 



(22) 



assuming a Gaussian and S ~ —k for simplicity. This 
estimate is valid for kvTo ^ where Tq is the initial 
temperature. 

For a fixed rj, all ensembles composed N < Nc par- 
ticles experience in a good approximation the same cool- 
ing time T ^ ojji/ri'^ for reaching the minimal tempera- 
ture (13). Nc is the critical particle number rendering 
the given rj critical. Note however, that this time scale 
is suboptimal for all ensemble sizes except N ^ Nf.. The 
scaling with TV for fixed laser power is more favourable as 
for standard cavity cooling [23]. Refer to fig. 6. For ini- 
tially stable ensembles the cooling rate is approximately 
the same until the instability point is reached. 

Let us briefly mention another aspect of the scaling of 
the Fokker-Planck equation (10) particularly useful for 
numerical simulations. Fixing \/N rj for different parti- 
cle numbers yields a cooling time r ^ TV/wr. Numerical 
simulations of the semiclassical equations (1) confirm this 
result, cf. figs. 7 and 8. There we have depicted the tem- 
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FIG. 5: (Colour on-line) Kinetic temperature for A'' = 500 and 
different pump strengths. Parameters: k — IOOljr and 5 = 
— K. The initial instability threshold is at 77 = lOr^c. Ensemble 
average over 5 initial conditions with 5 noise trajectories. 
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FIG. 6: (Colour on-line) Kinetic temperature for fixed pump 
strength and = {100, 500, 1000, 1750}, from bottom to top.. 
Parameters: 77 — 28ctjR, k = IOOojr and S = — k. The self- 
consistent threshold (14) is surpassed for all curves; the en- 
semble with N = 1750 is initially unstable. Ensemble average 
over 25 {N = 100) and 8 (higher particle numbers) initial 
conditions, respectively, and 5 noise trajectories. 



perature evolution for the pump strength being a fixed 
fraction of the critical value rjc for two different particle 
numbers, i. e. ^/Nirji = ^/7^r]2- In fig- 8 the threshold (7) 
is surpassed during the time evolution. 

In order to verify the Fokker-Planck equation (10), we 
consider the equation 



dt 



-A{F) 



dv. 



(23) 

for the kinetic temperature k^Tj^in ralv"^^ and close 
it by assuming {F{v,t)) to be a Gaussian velocity dis- 
tribution with temperature fceTkin for all times. This 
assumption is well justified for large detunings \5\ 3> 
and its predicted temperature reproduces the results ob- 
tained from the SDEs (1) quite well, cf. fig. 7. The main 
difference between the curves stems from the relatively 
small particle numbers in combination with pump val- 
ues close to threshold, where the hypothesis of separated 
time scales used to derive eq. (10) is no longer valid due 
to long-lived fluctuations. 
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FIG. 7: (Colour on-line) Comparison of the kinetic tempera- 
ture obtained from the SDEs (1) and as solution of eq. (23) 
for two different particle numbers and constant y/N-q. As 
expected, both ensembles evolve on a time scale ~ N/uj-r,. 
Parameters: NUo = -O.OIcjr, k = 100ci;r, Ac = NUo/2 - k, 
(i.e. 5 — —ft) and x/JV?? = SOoir = 0.8\/iVr;c. Ensemble av- 
erage of 50 (iV = 100) and 2 (iV = 1000) initial conditions, 
respectively, and 50 realisations of the white noise process. 
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FIG. 8: (Colour on-line) Temperature evolution above the 
self-consistent threshold (14). The homogeneous distribution 
is stable until the relation (7) is satisfied. This also limits the 
validity of eq. (23). Parameters: NUq = — a;R, k = IOOojr, 
Ac = NUo/2-K, (i.e. S = -k) and y/N-q = 200a;R = 2y^r]c. 
As \/7^ri = const., the photon number scales only ~ TV and 
not ~ A'^^ (superradiance effect). Furthermore, the equilib- 
rium temperature is the same for both ensembles. The sim- 
ulations were performed up to uj-Rt — 40A and revealed a 
temperature fcBlkin ~ 0.57k, which agrees very well with the 
theoretical prediction (16). Ensemble average of 50 {N = 250) 
and 25 {N = 1000) initial conditions, respectively, and 10 re- 
alisations of the white noise process. 



Of course, a thorough investigation of the validity of 
eq. (10) would require a numerical integration thereof. 
However, in steady state, all results of numerical sim- 
ulations of the SDE system (1) were found to be in 
excellent agreement with the predictions of the kinetic 
theory, both, below (e. g. fig. 2) and above (c. g. fig. 8) 
threshold. 



VIII. CONCLUSION AND OUTLOOK 

Collective light scattering from a dilute gas of cold par- 
ticles into a high-Q resonator mode under suitable con- 
ditions leads to friction forces and cooling of particle mo- 
tion even below the self-organisation threshold. In con- 
trast to standard cavity cooling the friction force below 
threshold only weakly depends on the particle number 
and leads to fast internal thermalisation towards a q- 
Gaussian velocity distribution of particles with average 
kinetic energy determined by the cavity linewidth only. 
Thus this constitutes a viable method for cooling very 
large ensembles of particles with sufficient polarisability 
independent of the need for a cyclic optical transition as 
well as a way to implement evaporative cooling for low 
densities, where no direct collisions occur. 

As only the polarisability and mass of a particle enter, 
one can easily envisage a combination of many species 
within the cavity to commonly interact with the same 
pump and cavity fields as a generalised form of sym- 
pathetic cooling without the need of direct interparticle 
interactions. The corresponding non-linear equation re- 
placing the Fokker-Planck equation (10) will then be of 
Balescu-Lenard type, cf. [21] and references therein. 
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